#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
#define ll long long
#define maxn 100+5
#pragma GCC optimize(2)
const double eps=1e-11;

void inputMartix(double a[][maxn],int n);
void printMartix(double a[][maxn],int N);
void LU(double a[][maxn],double b[],double L[][maxn],double U[][maxn],int N);
void solveY(double L[][maxn],double y[],double b[],int N);
void solveX(double U[][maxn],double x[],double y[],int N);
double Abs(double x);

int main(){
    double a[maxn][maxn];
    double L[maxn][maxn],U[maxn][maxn],b[maxn],x[maxn],y[maxn];
    int n;
    printf("请输入矩阵的阶数n,n<=100:");
    scanf("%d",&n);
    memset(L,0,sizeof(L));
    memset(U,0,sizeof(U));
    for(int i=1;i<=n;++i) L[i][i]=1;
    
    printf("请输入%d行%d列的矩阵:\n",n,n);
    inputMartix(a,n);//Ly=b,Ux=y;
    printf("请输入%d维的b向量:",n);
    for(int i=1;i<=n;++i)
        scanf("%lf",&b[i]);
    LU(a,b,L,U,n);
    
    printf("分解完成:\nL=\n");
    printMartix(L,n);
    printf("U=\n");
    printMartix(U,n);
    
    solveY(L,y,b,n);
    solveX(U,x,y,n);
    printf("由Ly=b,Ux=y,解得:中间向量y为:\n");
    for(int i=1;i<=n;++i)
        printf("y%d=%10.6lf\n",i,y[i]);
    printf("x解向量为:\n");
    for(int i=1;i<=n;++i)
        printf("x%d=%10.6lf\n",i,x[i]);
    return 0;
}

inline void LU(double a[][maxn],double b[],double L[][maxn],double U[][maxn],int N){
    int M[maxn];
    memset(M,0,sizeof(M));
    for(int k=1;k<=N;++k){
        double s[maxn];
        for(int i=k;i<=N;++i){
            double sum=0.0;
            for(int t=1;t<=k-1;++t)
                sum+=L[i][t]*U[t][k];
            s[i]=a[i][k]-sum;
        }
        double maxS=-1.0;
        int maxId;
        for(int i=k;i<=N;++i)
            if(Abs(s[i])>maxS){
                maxId=i;
                maxS=Abs(s[i]);
            }
        M[k]=maxId;
        if(maxId!=k){
            for(int t=1;t<=k-1;++t){
                double temp=L[k][t];
                L[k][t]=L[maxId][t];
                L[maxId][t]=temp;
            }
            for(int t=k+1;t<=N;++t){
                double temp=a[k][t];
                a[k][t]=a[maxId][t];
                a[maxId][t]=temp;
            }
            double temp=s[k];
            s[k]=s[maxId];
            s[maxId]=temp;
        }

        U[k][k]=s[k];
        if(k<N){
            for(int j=k+1;j<=N;++j){
                double sum=0.0;
                for(int t=1;t<=k-1;++t)
                    sum+=L[k][t]*U[t][j];
                U[k][j]=a[k][j]-sum;
            }
            for(int i=k+1;i<=N;++i){
                L[i][k]=s[i]/U[k][k];
            }
        }
    }
    for(int k=1;k<=N-1;++k){
        int t=M[k];
        double temp=b[k];
        b[k]=b[t];
        b[t]=temp;
    }
}

inline void solveY(double L[][maxn],double y[],double b[],int N){
    y[1]=b[1];
    for(int i=2;i<=N;++i){
        double sum=0.0;
        for(int j=1;j<=i-1;++j)
            sum+=L[i][j]*y[j];
        y[i]=b[i]-sum;
    }
}

void solveX(double U[][maxn],double x[],double y[],int N){
    x[N]=y[N]/U[N][N];
    for(int i=N-1;i>=1;--i){
        double temp=0.0;
        for(int j=i+1;j<=N;++j){
            temp+=U[i][j]*x[j];
        }
        temp=y[i]-temp;
        x[i]=temp/U[i][i];
    }
}


inline void inputMartix(double a[][maxn],int N){
    for(int i=1;i<=N;++i)
        for(int j=1;j<=N;++j)
            scanf("%lf",&a[i][j]);
}

inline void printMartix(double a[][maxn],int N){
    for(int i=1;i<=N;++i){
        for(int j=1;j<=N;++j)
            printf("%10.6lf ",a[i][j]);
        printf("\n");
    }
}



inline double Abs(double x){
    return x>0.0?x:-1.0*x;
}